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ABSTRACT 

We construct evolutionary models of the populations of active galactic nuclei (AGN) and supermassive black 
holes, in which the black hole mass function grows at the rate implied by the observed luminosity function, 
given assumptions about the radiative efficiency and the luminosity in Eddington units. We draw on a variety 
of recent X-ray and optical measurements to estimate the bolometric AGN luminosity function and compare to 
X-ray background data and the independent estimate of Hopkins et al. (2007) to assess remaining systematic 
uncertainties. The integrated AGN emissivity closely tracks the cosmic star formation history, suggesting that 
star formation and black hole growth are closely linked at all redshifts. We discuss observational uncertainties 
in the local black hole mass function, which remain substantial, with estimates of the integrated black hole 
mass density p, spanning the range 3 —5.5 x IO^'MqMpc^^. We find good agreement with estimates of 
the local mass function for a reference model where all active black holes have a fixed efficiency e = 0.065 
and Lboi/i'Edd ~ 0.4 (shifting to e = 0.09, Lboi/i'Edd ~ 0.9 for the Hopkins et al. luminosity function). In 
our reference model, the duty cycle of lO^M© black holes declines from 0.07 at z = 3 to 0.004 at z = 1 and 
10""* at z = 0. The decline is shallower for less massive black holes, a signature of "downsizing" evolution in 
which more massive black holes build their mass earlier. The predicted duty cycles and AGN clustering bias 
in this model are in reasonable accord with observational estimates. If the typical Eddington ratio declines at 
z < 2, then the "downsizing" of black hole growth is less pronounced. Models with reduced Eddington ratios 
at low redshift or black hole mass predict fewer low mass black holes (M, < 10** Mq) in the local universe, 
while models with black hole mergers predict more black holes at M, > lO^M©. Matching the integrated 
AGN emissivity to the local black hole mass density implies e = 0.075 x {p,/4.5 x 10^ M0 Mpc^-')^' for our 
standard luminosity function estimate, or 25% higher for Hopkins et al.'s estimate. It is difficult to reconcile 
current observations with a model in which most black holes have the high efficiencies e « 0. 16 — 0.20 predicted 
by MHD simulations of disk accretion. We provide electronic tabulations of our bolometric luminosity function 
and our reference model predictions for black hole mass functions and duty cycles as a function of redshift. 
Subject headings: cosmology: theory - black hole: formation - galaxies: evolution - quasars: general 



1. EVITRODUCTION 

The long-standing hypothesis that quasars are powered 
by accretion onto supermassive black holes (Salpeter 1964; 
Lynden-Bell 1969; Rees 1984) is now strongly supported by 
many lines of evidence, including the apparent ubiquity of 
remnant black holes in the spheroids of local galaxies (Rich- 
stone etal. 1998). The strong correlations between the masses 
of central black holes and the luminosities, dynamical masses, 
and velocity dispersions of their host spheroids (e.g., Magor- 
rian et al. 1998; FeiTai-ese & Memtt 2000; McLure & Dun- 
lop 2002; Marconi & Hunt 2003; Hiiring & Rix 2004; Fer- 
rarese & Ford 2005; Greene & Ho 2006; Graham 2007; Hop- 
kins et al. 2007b; Shankar & Ferrarese 2008; Shankar et 
al. 2008) imply that the processes of black hole growth and 
bulge formation are intimately linked. Theoretical models 
typically tie black hole growth to episodes of rapid star forma- 
tion, perhaps associated with galaxy mergers, and ascribe the 
black hole-bulge correlations to energy or momentum feed- 
back from the black hole or to regulation of black hole growth 
by the bulge potential (e.g.. Silk & Rees 1998; Kauffmann 
& Haehnelt 2000; Cavaliere & Vittorini 2000; Granato et al. 
2004, 2006; MuiTay, Quataert, & Thompson 2004; Cattaneo 
et al. 2005; Miralda-Escude & KoUmeier 2005; Monaco & 
Fontanot2005; Croton et al. 2006; Hopkins et al. 2006a, Mal- 
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bon et al. 2006). These correlations also make it possible to 
estimate the mass function of black holes in the local universe 
(e.g., Salucci et al. 1999; Yu & Tremaine 2002; Marconi et 
al. 2004; Shankar et al. 2004, S04 hereafter). This local mass 
function provides important constraints on the co-evolution of 
the quasar and black hole populations. The most general and 
most well known of these constraints is the link between the 
integrated emissivity of the quasar population, the integrated 
mass density of remnant black holes, and the average radia- 
tive efficiency of black hole accretion (Soltan 1982; Fabian & 
Iwasawa 1999; Elvis et al. 2002). 

In the paper we construct self-consistent models of the 
quasar population using a method that can be considered 
a "differential" generalization of Soltan's (1982) argument. 
Given assumed values of the radiative efficiency and the Ed- 
dington ratio L/LEdd, the observed luminosity function of 
quasars at a given redshift can be linked to the average growth 
rate of black holes of the corresponding mass, and these 
growth rates can be integrated forward in time to track the 
evolution of the black hole mass function. This modeling ap- 
proach has been developed and applied in a variety of forms 
by numerous authors, drawing on steadily improving obser- 
vational data (e.g. Cavaliere et al. 1973; Small & Blandford 
1992; Salucci et al. 1999; CavaUere & Vittorini 2000; Yu 
& Tremaine 2002; Steed & Weinberg 2003, hereafter SW03; 
Hosokawa 2004; Yu & Lu 2004; Marconi et al. 2004; S04; 
Merloni 2004; Vittorini, Shankar & Cavaliere 2005; Tamura 
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et al. 2006; Hopkins et al. 2007). The consensus of recent 
studies is that evolutionary models with radiative efficien- 
cies of roughly 10% and mildly sub-Eddington accretion rates 
yield a reasonable match to the observational estimates of the 
remnant mass function. However, uncertainties in the bolo- 
metric luminosity function of active galactic nuclei (AGN, 
a term we will use to describe both quasars and less lumi- 
nous systems powered by black hole accretion) and in the lo- 
cal black hole mass function remain an important source of 
uncertainty in these conclusions. 

We begin our investigation by constructing an estimate of 
the bolometric AGN luminosity function. Our estimate starts 
from the model of Ueda et al. (2003, hereafter U03), based 
on data from several X-ray surveys, but we adjust its parame- 
ters based on more recent optical and X-ray data that provide 
more complete coverage of luminosity and redshift. Compar- 
ison to the X-ray background and to the independent lumi- 
nosity function estimate of Hopkins, Richards & Hernquist 
(2007; HRH07 hereafter) gives an indication of the remaining 
uncertainties, associated mainly with bolometric corrections 
and the fraction of obscured sources. In agreement with Mar- 
coni et al. (2004) and Merloni (2004), we find similar trends 
in the evolution of AGN emissivity and the cosmic star for- 
mation rate, supporting models in which the growth of black 
holes occurs together with the formation of stars in their hosts 
(e.g., Sanders et al. 1988; Granato et al. 2004; Hopkins et 
al. 2006a). We also reassess current estimates of the local 
black hole mass function, finding that different choices and 
calibrations of the black hole-bulge correlation lead to sig- 
nificantly different results, with integrated mass densities that 
span nearly a factor of two. 

Our simplest models of the evolving black hole and AGN 
populations assume that all active black holes have a single 
radiative efficiency e and a single accretion rate m in Edding- 
ton units. The work of KoUmeier et al. (2006) suggests that a 
constant value of rh may be a reasonable approximation for lu- 
minous quasars, but the observational evidence on this point 
is mixed (e.g., McLure & Dunlop 2004; Vestergaard 2004; 
Babic et al. 2007; Bundy et al. 2007; Netzer & Trakhtenbrot 
2007; Netzer et al. 2007; Rovilos & Gorgantopoulos 2007; 
Shen et al. 2007b), and for lower luminosity, local AGN there 
is clearly a broad range of Eddington ratios (e.g., Heckmann 
et al. 2004). We also consider models in which m depends on 
redshift or black hole mass, and we consider a simple model 
of black hole mergers to assess their potential impact on the 
mass function. We will examine models with distributions of 
m values and a more realistic treatment of mergers in future 
work. 

Our modeling allows a consistency test of the basic sce- 
nario in which the observed luminosity of black hole accre- 
tion drives the growth of the underlying black hole population, 
and comparison to the local black hole mass function yields 
constraints on the average radiative efficiency and the typical 
accretion rate. For specified e and m, the model predicts the 
black hole mass function as a function of time, and combi- 
nation with the observed luminosity function yields the duty 
cycle as a function of redshift and black hole mass. These 
predictions can be tested against observations of AGN host 
galaxy properties and AGN clustering, and they can be used 
as inputs for further modeUng. While significant uncertain- 
ties remain, we find that a simple model of the black hole and 
AGN populations achieves a good match to a wide range of 
observational data. 



2. THE AGN BOLOMETRIC LUMINOSITY FUNCTION 

In order to constrain the accretion history of black holes, 
it is essential to know the shape and evolution with redshift 
of the AGN bolometric LF. The determination of this LF is 
not straightforward, as it must be made using the AGN LF 
observed in particular bands and converted into bolometric 
estimates using empirically calibrated bolometric corrections. 
The bolometric corrections have been shown by several au- 
thors to be redshift independent, with some possible trend 
with the intrinsic luminosity of the source (e.g., Elvis et al. 
1994; Marconi et al. 2004; Richards et al. 2006; Hopkins 
et al. 2007a). AGN surveys have been carried out at several 
energy ranges, and in each one the detection of sources can 
be highly affected by intrinsic or instrumental selection ef- 
fects. In different ranges of redshift and luminosity, the best 
statistics come from different wavebands or surveys, further 
complicating the effort to create a comprehensive AGN LF. 

An additional challenge in assembling a reliable and com- 
plete census of the AGN population is obscuration of the cen- 
tral engine by gas and dust, which may reside in the "molec- 
ular torus" of unified models (Antonucci 1993) or in the in- 
terstellar medium of the galactic host (Martinez-Sansigre et 
al. 2005; Rigby et al. 2006). Strong obscuration can elim- 
inate sources from surveys entirely, while uncorrected weak 
obscuration causes their luminosities to be underestimated. 
For solar composition gas, the absorption optical depth to X- 
ray photons is r = 2.QA{NH/\0^'^cm-~^){E /IksN)-^-'^ (e.g., 
Kembhavi & Narlikar 1999). For Galactic dust-to-gas ra- 
tio, the extinction in the rest-frame visual band is Ay = 
5.35mag x {Nh /IQ^^cmr-^) (e.g., Binney & Merrifield 2000). 
Hard X-ray selection is thus the least affected by obscura- 
tion, and deep X-ray surveys indeed reveal many faint AGNs 
that are missed by traditional optical selection criteria (e.g., 
Hasinger et al. 2005). Numerous studies suggest that the 
incidence of obscuration decreases towards high intrinsic lu- 
minosity, so that X-ray and optical LFs agree at the bright 
end but diverge at low luminosities (e.g., U03; Richards et al. 
2005; La Franca et al. 2005; Hasinger et al. 2005). Statis- 
tically speaking, there is fairly good correspondence between 
X-ray AGNs with logNn / cm < 22 and broad-line optical 
AGNs (e.g., Tozzi et al. 2006). Even 2-10 keV selection 
becomes highly incomplete for Compton-thick sources, with 
logNn I cm~^ > 24, and the best constraints on this population 
come from the normalization and spectral shape of the X-ray 
background. 

Given these considerations, we have chosen to base our 
bolometric LF estimate principally on the work of U03, who 
compiled a vast sample from Chandra, ASCA, and HEAO-l 
surveys and inferred the absorption-corrected LF in the rest- 
frame 2-10 keV band out to z ^ 3. We adopt the luminosity- 
dependent bolometric correction of Marconi et al. (2004); we 
include a dispersion of 0.2 dex in this correction (see, e.g., 
Elvis et al. 1994), but omitting the scatter would not signif- 
icantly affect our results. U03 fit their data with a parame- 
terized model of luminosity-dependent density evolution. We 
adopt the model and adjust its parameters in some redshift 
ranges to better fit other data sets, especially at high redshifts 
and at bright luminosities at low redshifts. We note that many 
other studies have reported general trends of the LF evolu- 
tion and obscuring columns similar to those of U03 (e.g.. La 
Franca et al. 2005), and that infrared surveys also suggest a 
substantial population of obscured AGNs at low luminosities 
(Treister et al. 2006; see also Ballantyne and Papovich 2007). 
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The number of sources per unit volume per dex of lumi- 
nosity logLx = logLi-iokeV is fitted by U03 using a double 
power-law multiplied by an "evolution" term: 



$(Lx,z) =e{z,Lx 



(1) 



where 



e[Lx,z)-< (i+,jP,t(i+2)/(i+,^(i^))]P. if z>Zc{Lx) 



with 



Zc{Lx) = 



0.335 



ifLx>La 
if Lx < La 



(2) 
(3) 



We tune the value of parameters in order to provide a good 
fit to the overall set of data presented in Figure [T] The 
full list of parameters is given in Table [T] Some of the pa- 
rameters assume different values in different redshift inter- 
vals, in which case we apply a linear interpolation in z be- 
tween these intervals, as reported in Table [1] The X-ray 
LF of equation ([TJ is converted into a bolometric LF using 
the fit to the bolometric correction given by Marconi et al. 
(2004), \ogL/Lx = 1. 54 + 0.24C + 0.012(2 -O.OOISC^ ^j^j^ 
C = logL/L0 — 12, and L© = 4 x lO-'-' ergs~', then con- 
volved with a Gaussian scatter of dispersion 0.2 dex to ac- 
count for dispersion in the bolometric correction. We find the 
bright-end slope 72 should vary with time to match the data, 
steepening at very low redshifts (z < 0.8) and at high (z > 4) 
redshifts. We caution that the "evolution" term e{Lx,z) sub- 
stantially modifies the shape of the luminosity function below 
L w lO'^^ergs^'. For example, the effective LF slope in the 
range 10'*^ergs^' < L < 10'*"^'^ergs^' changes from —1.4 at 
z ~ 2 to —0.7 at z even though we keep 71 (appearing in 
equation IT]) fixed at 0.86. 

Figure [T] compares our model of the AGN LF to a large 
collection of data from optical surveys (Pei 1995; Wisotzki 
1999; Fan et al. 2001, 2004; Kennefick et al. 1994; Hunt et 
al. 2004; Wolf et al. 2003; Richai'ds et al. 2005, 2006; Jiang et 
al. 2006; Cool et al. 2006; Bongiorno et al. 2007; Fontanot et 
al. 2007; Shankar & Mathur 2007) and X-ray surveys (Barger 
et al. 2003; Ueda et al. 2003; Bai-ger & Cowie 2005; Barger 
et al. 2005; La Franca et al. 2005; Nandra et al. 2005; Silver- 
man et al. 2008). Optical data have been converted from the 
Z?-band into bolometric quantities using L^^vbCb = L, where 
Li,^ is the monochromatic luminosity at vb = 6.8 x lO'^Hz 
(no 4400 A). We use an average value of Cb — 10.4 consis- 
tent with Richards et al. (2006). Note that Marconi et al. 
(2004) and Hopkins et al. (2007; HRH07 hereafter) suggest 
values 30%-50% lower depending on luminosity, while Elvis 
et al. (1994) proposed a higher value of 11.8. In Figure [1] 
we plot L$(L), instead of simply ^(L) itself, to highlight the 
luminosity bins where most of the energy is emitted. For the 
same reason, we will later plot the black hole mass function 
as M,<i>(M,) to highlight the mass bins in which most of the 
density resides. 

The dot-dashed curves in Figure [T] show our model LF in- 
cluding only sources with logA^/z/cm^^ < 22, which we cal- 
culate using the luminosity-dependent column density distri- 
bution functions of U03. The dashed curves show the contri- 
bution of all sources with logA^/^/cm^^ < 24, while the solid 
curves include higher column density (Compton-thick) sys- 
tems. We assume that the number of Compton-thick sources 
in each luminosity bin is equal to the number of sources in 



the column density range 23 < logA^/z/cm^^ < 24. We will 
henceforth follow the X-ray convention of describing systems 
with logA^n/cm^^ < 22 as "Type I" AGN and systems with 
22 < logNn /cm"2 < 24 as Type II AGN. We roughly expect 
the "Type I" curve (dot-dashed) to agree with optical survey 
data points and the "Type In-Type 11" curve to agree with X- 
ray survey data points. The agreement between the model and 
the data is overall fairly good, though in some regimes differ- 
ent data sets give seemingly incompatible results, making it 
difficult to judge the validity of the model itself. 

HRH07 have used their observationally constrained theo- 
retical models to estimate the fraction of AGN that are "miss- 
ing" from observational samples because of obscuration, as 
a function of waveband and luminosity. In Figure [21 we plot 
data points corrected for these obscured fractions, with the 
same model curves as Figure [H If the data, the obscuration 
corrections, and our model were all perfect, then all of the data 
points should line up with each other and with the solid model 
curves. The obscuration corrections do reduce the discrepan- 
cies among data sets, most notably in the range z ~ 0.7 — 3.2, 
but they do not remove all of the differences. The upper enve- 
lope of the data points is generally close to our model near the 
peak of L$(L), though at the highest luminosities the model 
exceeds the data (principally SDSS) for z ~ 0.7 — 2.0. 

In Figure [3 we show the integrated intensity obtained from 
our model AGN LF compared with all the available data on 
the cosmic X-ray Background (XRBG) for energies above 1 
ke V. The data are a collection of old and new results presented 
by Frontera et al. (2007), plus the recent results from INTE- 
GRAL (Churazov et al. 2006). While the shape of the XRBG 
is now well established at low energies (2-10 keV band), dif- 
ferent studies imply normalizations that differ by up to 40%. 
The minimum estimate for the normalization was found by 
the very first HEAO experiments (e.g., Marshall et al. 1980; 
filled circles in Figure [3]l, while the much more recent Chan- 
dra and XMM results point towards higher values (shaded 
bands in Figure [3l see Comastri 2004 for a review). Sev- 
eral groups have assumed that the HEAO measurements were 
correct in shape but underestimated in normalization because 
of flux calibration and/or instrumental background subtrac- 
tion (Frontera et al. 2007), implying a true XRBG spectrum 
that is 1.4 times the HEAO spectrum at all energies. How- 
ever, recent measurements from INTEGRAL (Churazov et al. 
2006; open circles in Figure [3]l and from the PDS instrument 
aboard BeppoSAX (Frontera et al. 2007) produce results close 
to the original HEAO measurements in the ~ 20-100 keV 
range. Based on the agreement of multiple experiments in 
overlapping energy ranges, we tentatively conclude that the 
true XRBG spectrum approximately follows the INTEGRAL 
points over the range ^ 5-100 keV. 

The solid curve in Figure [3 shows the XRBG predicted 
by our AGN bolometric LF model with the U03 column 
density distribution. Following U03, we use the PEXRAV 
code (Magdziarz & Zdziarski 1995) to generate spectra of 
families of AGN with different column densities, also tak- 
ing into account Compton down-scattering for highly ob- 
scured sources (Wilman & Fabian 1999). We compute 
the spectra of Compton-thick (logA^/^/cm^^ > 24) sources 
assuming logA^/z/cm^^ — 24.5. Dashed, dot-dashed, and 
triple-dot-dashed curves show the cumulative contributions 
of sources with logA^///cm^- < 22,23, and 24, respectively. 
The Compton-thick sources increase the predicted XRBG by 
a factor of about 1.3 at £ > 20 keV, with smaller contribu- 
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tions at lower energies and at much higher energies. Our full 
model would be a good match to the HEAO spectrum renor- 
malized by a constant factor of ^ 1.4 as proposed in some 
earlier studies. However, it significantly overpredicts the IN- 
TEGRAL and SAX/PDS results at £ > 10 keV; these ai-e well 
matched by the contribution of logA^///cm^^ < 24 sources 
alone, with no Compton-thick contribution. As noted above, 
we employ U03's prescription for the luminosity dependence 
of the Compton-thick fraction. If we use a simpler model 
in which the number of Compton-thick sources is 40% of the 
number of log A^n/cm^^ < 24 sources at all luminosities, then 
we obtain very similar XRBG predictions (higher by a few 
percent near the peak and even closer at other energies). 

We conclude from Figure[3]that the Compton-thick fraction 
in our standard model is probably an upper limit on the true 
fraction, and in our calculations below we will also consider 
a model in which the Compton-thick fraction is zero. How- 
ever, Gilli et al. (2007) have proposed a successful model to 
fit the INTEGRAL measurements of the XRBG that uses a 
combined contribution from unobscured sources (which they 
define to be only those with logA^///cm^^ < 21) a factor 2 
below our estimate and a higher contribution from obscured 
sources up to logA^///cm^^ —26, which they derive by assum- 
ing that these sources have a Gaussian distribution of spectral 
indices around a mean slope (F) = 1.9. With this result in 
mind, we will also explore models with a higher fraction of 
Compton-thick sources (but our standard estimate of unob- 
scured sources). 

Reproducing the high background normalization found by 
XMM in the 2-5 keV range would require a substantially 
higher contribution from Type I AGN, since obscured AGN 
have little emission at these energies. Such an increase seems 
implausible on other grounds. Following Schirber & Bullock 
(2003) and Miralda-Escude (2003), we have estimated the 
average hydrogen ionization rate produced by Type I X-ray 
AGN with optical luminosity Mb < — 15 assuming UV spec- 
tral energy index auv = 1.57. Using the redshift-dependent 
X-ray-to-optical ratio by Steffen et al. (2006), we find that 
the Type I AGNs produce an emissivity that saturates the es- 
timated UV background intensity at z < 1 . We therefore con- 
clude that the true AGN contribution to the XRBG is probably 
closer to our model prediction than to the XMM band shown 
in Figure [3] If future missions confirm a high normalization 
of the XRBG in the 2-5 keV range, it would probably mean 
that a substantial contribution from X-ray binaries in normal 
galaxies is present at these energies. On the other hand, we 
note from Figure [T] that our AGN LF defined for sources with 
logA'/z/cm^ < 22 is consistent with, or even lower than, all 
optical survey estimates, so we do not expect a much lower 
contribution from these sources. 

HRH07 have recently undertaken an exercise similar to 
ours, attempting to construct a bolometric AGN LF that is 
consistent with a wide range of data over a large span of 
energies and wavelengths. Figure |4] compares our model to 
theirs. We find good agreement in the shape and overall nor- 
malization of the AGN LF in the redshift range z < 1 and 
3.5 < z < 4.5, while at intermediate redshifts our model LF 
is systematically lower The normalization difference arises 
because HRH07 adopt an X-ray bolometric correction that is 
about 30% higher (at typical luminosities) than the Marconi 
et al. (2004) bolometric correction used here. If we adopt 
the HRH07 bolometric correction, then we find good agree- 
ment in the range 1 < z < 3, but our model LF is higher at 
high and low redshifts. At z > 4 we adopt a steeper bright- 



end LF slope based on the results of GOODs (Fontanot et al. 
2006), Cool et al. (2006) and Shankar & Mathur (2007), who 
find evidence for a steeper slope relative to the estimates of 
Richards et al. (2006). However, at 1 < z < 2.5 the HRH07 
bright-end slope is steeper than ours, since they effectively re- 
quire a match to the SDSS-based measurements of Richards 
et al. (2004) in this redshift range, while we retain the shal- 
lower slope adopted by U03, also supported by recent X-ray 
LF determination of Silverman et al. (2008). Recent Spitzer 
studies (Hickox et al. 2007; Polletta et al. 2008) suggest a sig- 
nificant fraction of obscured AGN even at high luminosities, 
leaving some room for a slope difference between the opti- 
cal and bolometric LFs. We regard the differences between 
our model LF and that of HRH07, derived from independent 
attempts to match the full range of current data, as a reason- 
able representation of the remaining systematic uncertainties 
in determination of the bolometric AGN luminosity function. 
We will show in § 14.31 that the difference between these two 
determinations does not alter our overall conclusions, but the 
differences in normalization and shape at intermediate red- 
shifts do lead to different preferred accretion parameters for 
the black hole population. 

3. THE LOCAL BLACK HOLE MASS FUNCTION AND THE 
INTEGRATED SOLTAN ARGUMENT 

3.1. The local black hole mass function 

Marconi et al. (2004) and S04 found similar results for 
the local black hole mass function, using somewhat different 
methods. Both groups found that starting from the relation be- 
tween black hole mass and bulge velocity dispersion (M, — a) 
or the relation between black hole mass and bulge luminosity 
(M, — Lsph) led to similar local mass functions. However, a 
more recent analysis by Lauer et al. (2007a) and Tundo et 
al. (2007) argues that these two methods do not provide the 
same answer We therefore revisit some of the uncertainties 
in computing the local mass function. 

Figure |5] compares the local mass functions obtained using 
different relations between black hole mass and host galaxy 
properties. Following S04, we start with the galaxy LF in 
the r* band by Nakamura et al. (2002), who used light con- 
centration parameters to distinguish early and late type galax- 
ies and derived LFs for both. After correcting the Petrosian 
magnitudes to total magnitudes by adding —0.2 mag, we con- 
vert the LF into a LF of spheroids, which we compute as fol- 
lows. First we compute the numerical fraction of EUiptical-SO 
(spirals Sab-Sbc-Scd) in the early-type (late-type) galaxy LF, 
then correct each morphological galaxy type for its respec- 
tive spheroidal luminosity component, as given in Table 1 of 
Fukugita et al. (1998) for the r-band (see also Yu & Tremaine 
2002 and Marconi et al. 2004). Note that this method is equiv- 
alent to correcting the luminosities themselves adopting an 
average weighted fraction of /sph = 0.83 — 0.85 for early-type 
galaxies, and /sph = 0.27 — 0.3 for late-type galaxies, as done 
in S04 (and references therein). 

The solid line in Figure |5] shows our result using the 
M, — Lsph relation calibrated by McLure & Dunlop (2002; 
their equation 6) for a sample of 18 black holes in inactive 
galaxies, where we correct the magnitudes for our cosmol- 
ogy. Using the calibration in equation (A9) of Tundo et al. 
(2006) would give a similar result. These calculations yield 
a black hole mass function ^ 34% higher than that in S04 
(shown with solid squares). We include a Gaussian scatter of 
0.3 dex around the mean M, — Lsph relation in both cases. The 
intrinsic scatter to insert in the local relations between black 
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hole mass and host galaxy properties is empirically uncertain 
because it is similar in magnitude to the observational uncer- 
tainties themselves. However, most authors estimate a scatter 
of about 0.3 dex, which is close to the value predicted in the 
numerical simulations of Hopkins et al. (2007b). In the fol- 
lowing we will always adopt this value of the scatter for all 
our computations of the local mass function. 

If we instead use the galaxy velocity dispersion function 
and the M, — a relation in equation (A5) of Tundo et al. 

(2006) , with a scatter of 0.3 dex, we get the short-dashed line 
in Figure |5] close to the central estimate of S04. Here we use 
the velocity dispersion function by Sheth et al. (2003), which 
includes their estimate for the contribution of the bulges in 
spirals, which in turn accounts for ^25% of the total esti- 
mated black hole mass density. The Sheth et al. (2003) esti- 
mate of the velocity dispersion function is in very good agree- 
ment with the velocity function estimated by S04 through the 
bivariate relation of galaxy luminosity and velocity disper- 
sion. However, using the M, — a relation in Marconi et al. 
(2004) yields the dot-dashed line, which is higher than the 
S04 determination by about 1 — cr. The M, — a relation given 
in Ferrarese & Ford (2005) is steeper than the Tundo et al. 

(2007) relation, and adopting it yields a similar mass den- 
sity to S04, but a local mass function shifted towards higher 
masses, shown with a triple-dot dashed line in Figure |5] (see 
also Wyithe 2004). Both of these estimates yield a local mass 
function below the M, — Lsph-based local mass function esti- 
mate at M, < 10^'^ M0, as also noted by S04. 

In the following we will adopt the grey band of Figure |5] 
which spans the range of these estimates, as representative 
of the mean and the systematic uncertainties of present esti- 
mates of the local mass function. The integrated mass den- 
sity of the local black hole population is p, = (3.2 — 5.4) x 
lO^MoMpc^^ (for h — 0.7). Figure |5] also presents three 
additional estimates of the local mass function. Stars show 
the results of combining the Bell et al. (2003) galaxy bary- 
onic mass function (corrected for spheroid fraction as above) 
with the Haring & Rix (2004) estimate of the relation be- 
tween black hole mass and spheroid stellar mass, again as- 
suming 0.3 dex intrinsic scatter This result is in reasonable 
agreement with the M, — a estimates, though uncertainties in 
the stellar mass-to-light ratio are a remaining source of sys- 
tematic uncertainty. The dotted curve shows the estimate of 
Hopkins et al. (2007b), based on a "fundamental plane" re- 
lation between black hole mass, galaxy velocity dispersion, 
and spheroid half-light radius; it is more sharply peaked than 
the grey band and implies a higher integrated black hole mass 
density. Open circles show Graham et al.'s (2007) estimate 
based on the correlation between black hole mass and the Ser- 
sic light-profile index of the galaxy spheroid. This estimate is 
similar to Hopkins et al.'s (2007b) near the peak, but it im- 
plies an even sharper fall off towards low black hole masses. 
Further discussion of the discrepancies in local mass function 
estimates, and possible routes to alleviate them, will appear in 
Shankar & Ferrarese (in preparation). 

3.2. The integrated mass density 

Before turning to the evolution of the differential black hole 
mass function, the central theme of this paper, we briefly re- 
visit the classic Soltan (1982) argument, which relates the 
integrated black hole density to the integrated emissivity of 
the AGN population. If the average efficiency of converting 
accreted mass into bolometric luminosity is e = L/MinflowC^, 



where MinOow is the mass accretion rate, then the actual accre- 
tion onto the central black hole is M, = (1 — e)Mi„flow, where 
the factor 1 — e accounts for the fraction of the incoming mass 
that is radiated away instead of being added to the black hole. 
The rate at which mass is added to the black hole mass func- 
tion is then given by 

^=^—^rmLd\ogL. (4) 
dt ec^ Jq 

The mass growth rate implied by equation ^ and our esti- 
mate of the AGN LF from §2 is shown by the solid line in Fig- 
ure|6h. We set the radiative efficiency to a value of e = 0.075, 
as it provides a cumulative mass density in agreement wit h the 
median estimate of the local mass density discussed in § 13.11 
At each time step we integrate equation (|4|i down to the ob- 
served faint-end cut in the 2 — 10 keV AGN LF, which we 
parameterize as 

l0gLMIN,2-I()kev(z) = logLfl, 2- lOkeV + 2.5 log( 1 + z) ■ (5) 

We set logL2-i{)kev/(ergs^') ~ 41.5, in agreement with the 
faintest low redshift objects observed by U03 and La Franca 
et al. (2005). For a typical Lopt/Lx, equation ^ yields an 
optical luminosity of Mg ^ —22 at z ^ 6, comparable to the 
faintest AGN sources observed by Barger et al. (2003) in 
the 2 Msec Chandra Deep Field North (see also Figure [T] and 
Shankar & Mathur 2007). At each time step we compute the 
minimum observed luminosity given in equation (|5]l and con- 
vert it into a bolometric quantity Lmin (z) applying the adopted 
bolometric correction by Marconi et al. (2004). 

Dashed and dot-dashed lines in Figure [6^ show two recent 
estimates of the cosmic star formation rate (SFR) as a function 
of redshift, from Hopkins & Beacom (2006), reported with its 
3-(T uncertainty region (dark area), and Fardal et al. (2007). 
We have multiplied both estimates by a redshift-independent 
factor of 8 X 10^^. Since local estimates imply a typical ra- 
tio M,/Mstar 1.6 X 10^^ for spheroids (e.g., Haring & Rix 
2004), this is a reasonable scaling factor if roughly 50% of star 
formation goes into spheroidal components (see also Marconi 
et al. 2004 and Merloni et al. 2004). The agreement between 
the inferred histories of black hole growth and star formation 
suggests that the two processes are intimately linked. In par- 
ticular this association seems to hold down to the last sev- 
eral Gyrs, even at z < 1 .5 when disk galaxies are expected to 
dominate the SFR. A possible link between black hole growth 
and star formation in disks could arise from re-activations in- 
duced by tidal interactions between satellite and central galax- 
ies (e.g., Vittorini et al. 2005). Also, bars could possibly fun- 
nel gas into the central black hole, though empirical studies 
cast some doubt on this mechanism as a primary trigger for 
black hole growth (Peeples & Martini 2006 and references 
therein). 

Figure IS]? presents the same comparison in integrated form 
(see also De Zotti et al. 2006 and Hopkins et al. 2006b). Solid 
squares show the black hole mass density obtained by con- 
verting the z— I and z — 2 galaxy stellar mass function into a 
black hole mass function by assuming a ratio M,/Mstar equal 
to the local one (i.e., 1.6 x 10^^). The galaxy stellar mass 
function has been computed from the Caputi et al. (2006) K- 
band galaxy luminosity function, assuming an average mass- 
to-light ratio Ms,ar/LA:=0.4 at z = 1 and Mstai/i;f=0.3 at z = 2. 
The latter values have been obtained from the Pegase2 code 
(Fioc & Rocca-Volmerange 1997) by taking a short burst of 
star formation (< 10^ yr) and a Kennicutt double power-law 
stellar Initial Mass Function. The quoted values for M^ii^/Lk 
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can be taken as lower limits, as other choices of the parame- 
ters in the code would tend to increase their value. However, 
we also note that our result on the stellar mass function is 
in good agreement with the recent estimate by Fontana et al. 
(2006). Our scaling factor of 1.6 x lO^"' implicitly assumes 
that all the stellar mass in the luminous galaxies probed by 
these high-redshift observations resides in spheroidal compo- 
nents today, and is therefore associated with black hole mass. 

Figure|6]suggests that the ratio of black hole growth to SFR 
is approximately the same at all redshifts, and suggests a close 
link between black hole accretion and star formation. If the 
average ratio of black hole mass to stellar mass were much 
higher than the local value at z = 1—2, then the squares 
in Figure |6}3 would shift above the black hole mass density 
curve p,(z). Increasing M^tm/LK would move the squares 
higher still, while associating only a fraction of the high- 
redshift stellar mass with present day spheroids would move 
them lower Our conclusions are in marginal agreement with 
those of McLure et al. (2006) and Shields et al. (2006) and in 
some disagreement with those of Peng et al. (2006), who find 
M,/Mstar in gravitationally lensed quasar hosts at z 1 — 2 a 
factor ^ 3 above the local value. However, several observa- 
tional biases may effect studies of the M,/Mstar ratio in high 
redshift quasar samples (Lauer et al. 2007b). Uncertainties in 
the local normalization of p, and in the stellar mass-to-light 
ratios still leave a fair amount of wiggle room, but our results 
in Figure|6]disfavor models in which black holes "grow first" 
and spheroid star formation "catches up" at later times. Al- 
ternatively, if the average efficiency e of black hole accretion 
is lower at high redshifts, then the p, {z) curve in Figure 
would shift upwards at these redshifts accordingly. 

4. EVOLVING THE BLACK HOLE MASS FUNCTION 

4.1. Method 

Our goal is to calculate the evolution of the black hole mass 
function implied by the bolometric AGN LF described in §2, 
given assumptions about the radiative efficiency and the typ- 
ical accretion rate. In the following we will use the symbol 
<l>(x) to denote mass and luminosity functions in logarithmic 
units of L or M,, i.e. 



^{x)=n{x)x\\\{\Q), 



(6) 



where n(x) is the comoving space density of black holes in 
the mass or luminosity range x — > .jc + dx, in units of Mpc^-' 
for h = 0.7. 
We define the Eddington accretion rate to be 



Edd ■ 



^Edd 
0.1 c2 



:22 



f M, 



V 10'' Mc 



M, 



©yr 



(7) 



and the dimensionless accretion rate m = M./Mndd, where 
LEdd is the standard Eddington luminosity (for Thomson scat- 
tering opacity and pure hydrogen composition) at mass M,. 
The black hole growth rate M, is related to the large scale 
accretion rate by M, = (1 — e)Minflow, where L ~ eMingowC^, 
because a fraction e of the mass is radiated away before enter- 
ing the black hole. We define / = e/(l — e), and /o.i = //O.l, 
so that a black hole of mass M, growing at a dimensionless 
rate m has bolometric luminosity 



L = eMinflowC^ = 0.1/o.iM.c- = fo.imlM, , 



(8) 



where / = LEdd/M, = 1.26 x 10'"*ergs"' Mq . Note that we 
define Mndd for / = 0.1, so the mass Eddington ratio m is 



Hnked to the luminosity Eddington ratio A by the radiative ef- 
ficiency, i.e. 

A=-^=m/o.i. (9) 

i-Edd 

A black hole accreting at a constant value of m grows expo- 
nentially in time with a timescale fg,owth — M,/M, = tg/m = 
4.5 X lO^m^'yr, where is equal to the Salpeter (1964) 
timescale for /o.l — 1. 

The evolution of the black hole mass function n{M,,t) is 
governed by a continuity equation (Cavahere et al. 1973, 
Small & Blandford 1992) 



dn 
'di 



d{M,{m)n{M„t)) 



(10) 



where (m) is the mean dimensionless accretion rate (averaged 
over the active and inactive populations) of the black holes of 
mass M, at time t. This evolution is equivalent to the case in 
which every black hole grows constantly at the mean accre- 
tion rate (m). In practice, individual black holes turn on and 
off, and there may be a dispersion in m values, but the mass 
function evolution depends only on the mean accretion rate as 
a function of mass. 

All models in this paper assume a single accretion rate 
m — niQ. At any given time, a black hole is either accreting 
at iriQ or not accreting. In some models we allow the char- 
acteristic accretion rate niQ to depend on z, or on M,. The 
assumption of a single lii is clearly not valid for low luminos- 
ity AGNs in the nearby universe, which have a wide range of 
Eddington ratios (e.g., Heckman et al. 2004, Greene & Ho 
2007). However, Kollmeier et al. (2005) find that luminous 
AGN at 0.5 < z < 3.5 have a narrow range of Eddington ra- 
tios, with a peak at A ^ 1 /4 and a dispersion of 0.3 dex. Since 
this dispersion includes contributions from random errors in 
black hole mass estimates and bolometric corrections, the true 
dispersion should be even smaller. Netzer et al. (2007) find a 
similar result, with a slightly larger dispersion, from a sample 
centered at z 2.5. We will consider models with multiple m 
values in future work, but single-;7t models (also adopted by, 
e.g., Marconi et al. 2004 and S04) are a good starting point 
for understanding black hole growth. 

If there is a single accretion rate mo, then the duty cycle of 
black hole activity (i.e., the probability that a black hole of 
mass M, is active at a particular time) is given by the ratio of 
luminosity and mass functions. 



PoiM„z) = 



HL,z) 



i/logM. 



<1, 



(11) 



where M, — L/(/o.imo/) is the black hole mass that corre- 
sponds to luminosity L. Models with constant niQ and e have 
L oc M,, making the Jacobian factor unity. A physically con- 
sistent model must have Pq {M, , z) < 1 ; there must be enough 
black holes to produce the observed luminous AGNs. 

Our strategy is to start with an assumed black hole mass 
function n(M,,z,) at an initial redshift z/, then track the char- 
acteristic curves M,(M,., ,z) of equation ( fTOb by direct inte- 
gration 

r dt 
M,{M,,i,z)= / {m,)M,{z')—dz' 



moMM.,z')M,{z')^dz' . 

dz 



(12) 
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Here Po{M,,z') is given by equation (fTTT l with the observed 
luminosity function, and the evolved n{M,,z) is given by 

n{M„z)^niM.j,Zi)^^, (13) 

where M, is the black hole mass that corresponds to initial 
mass M, i. To put this calculation in more physical terms: we 
take advantage of the fact that equation ( fTOl i is equivalent to 
having every black hole grow at the rate {m{M,,z))- Starting 
with a set of logarithmically spaced initial values of M, ,, we 
integrate the masses forward in time with a mid-point scheme 
from redshift step Zj to step Zj+] = Zj — Az 

M,j+i = M,j + moPo{M,j+i/2,Zj+i/2)M,j+i/2^^Az , 

(14) 

where the values at the mid-point j + 1/2 are evaluated by 
extrapolating black hole masses for a half step Az/2 using 
dM,/dt at the beginning of the step and the duty cycle Pq is 
computed using the mid-step mass function and the luminos- 
ity function evaluated at Zy+1/2 = zj + Az/2. We choose suffi- 
ciently small redshift steps Az such that the total mass density 
added at each iteration matches the one obtained by integra- 
tion of the bolometric luminosity function (equation |4]i, and 
we have confirmed that smaller Az yields essentially iden- 
tical results. Since the number of black holes is conserved, 
n (M. , z) AM. =n (M.,,- , z) AM.,,- . 

As an aside, we note that this approach is mathematically 
equivalent to solving the continuity equation in the form 

dn{M„t) ^ mo d<i>{L,z) 

dt f,ln(10)2M. dlogL ^ ' 

used by, e.g. Small & Blandford (1992) and Marconi et al. 
(2004). Equation ( fTSl ) follows from equation (fTOb with the as- 
sumption of a single mo and /o 1 , and it is also valid in the case 
of a redshift dependent mo. From equation ( fTsT i, it is explic- 
itly clear that the evolution of the black hole mass function 
depends only on the initial conditions, the observed luminos- 
ity function, and the values of mo and fo.\ =L/ (M./mo), given 
the assumptions made here. The Soltan argument (§ 13. 2b re- 
lates the integrated black hole mass density to the integrated 
emissivity of the AGN population. The continuity equation 
yields the full evolution of the black hole mass function in 
terms of the full luminosity function; in essence, it applies the 
Soltan argument as a function of mass. We have checked that 
our method yields consistent results against direct integration 
of the equation. 

For our integrations, we generally start at z, = 6 and de- 
termine n(M. ,,z,) from equation ( fl4l l with our estimate of 
the luminosity function at z = 6 and an assumed Po ~ 0.5. 
This relatively high initial duty cycle implies we start with 
nearly the minimal black hole mass function required to re- 
produce the luminosity function. By z ~ 3.5 the integration 
has essentially forgotten the initial value of Pq unless we set 
it much lower (e.g., Po < 0.01), since the accreted mass is 
much larger than the mass in the initial black hole population. 
The clustering of luminous AGN at z 3 — 4 implies duty 
cycles of at least several percent (Shen et al. 2007; Shankar 
et al., in preparation). In our standard calculations, we set Po 
to zero for masses M..min(z) that would produce AGN lumi- 
nosities below the minimum observed luminosity at that red- 
shift (indicated by the cutoffs in Figure[T]i; we do not extrap- 
olate the luminosity function below these minimum values. 
Since Lmin (z) drops towards lower z, we must have a supply 



of black holes in place to provide AGNs of lower luminos- 
ity as these become visible. To ensure this, we set our initial 
black hole mass function to 

n{M,^i,Zi) = 8 X n{M,MiN{zi),Zi) x — — 

(16) 

i.e., below M. minI-Z/) we boost n{M,,Zi) by a factor of 8 
and add an M^^^ rise (see Figure |2l below). This high ini- 
tial abundance of low mass black holes ensures that we are 
never forced to duty cycles Po > 1 . Once M. significantly ex- 
ceeds M, min(z), the black hole mass function is dominated 
by growth rather than initial values, and «(M,,z) is insensi- 
tive to the assumed «(M, ,,z,). We adopted this procedure so 
that our results would not depend on the extrapolation of the 
luminosity function into regions where it is not observed. In 
practice, we find that extrapolating the luminosity function as 
a power law down to very faint luminosities yields similar re- 
sul ts. W e f ollow the Lmin procedure for our standard m odels in 
§§|42]and |431 and use extrapolation of the LF in §§|43E6] 
where it gives more stable numerical solutions. 

4.2. The reference model 

Figure |7] shows the results of our calculations for a refer- 
ence model with our standard estimate of the AGN bolomet- 
ric luminosity function and accretion parameters m = 0.60 
and e = 0.065 (/o.i ~ 0.7), which yield good overall agree- 
ment with the average estimate of the local black hole mass 
function. Panel |7^ shows the evolution of the mass function 
starting from the initial condition at z = 6 (solid triangles) 
and continuing through to z = 0. Because of the luminosity- 
dependent density evolution in the observed luminosity func- 
tion, the massive end of the black hole mass function builds up 
early, and the lower mass regime grows at later redshifts. For 
M, > 10^^ M0, the mass function is almost fully in place by 
z = 1 . Panel|7j5 plots M. $ (M, ) , proportional to the fraction of 
black hole mass per logarithmic interval of M., which allows 
better visual comparison to the observed local mass function 
and highlights the contribution of each black hole mass bin 
to the total mass density at each time. The accretion model 
agrees well with our estimate of the local mass function (grey 
band) up to 10^'^ Mq. At higher masses our model exceeds 
the observational estimate, but in this regime the estimate re- 
lies on extrapolation of the scaling relations, and it is sensitive 
to the assumed intrinsic scatter In addition, the high-mass end 
of the predicted mass function is sensitive to the bright end of 
the LF at z < 2 (see § l4.3l l. The open circles with error bars 
show the estimate of the local mass function by Graham et al. 
(2007), which we cannot reproduce even approximately with 
constant mo- If this estimate is correct, then the low end of the 
luminosity function must be produced mainly by high mass 
black holes accreting at low m, so that the predicted growth of 
low mass black holes is reduced (see § l4.5l l. 

Figure|7}; plots the duty cycle as a function of mass for dif- 
ferent redshifts, as labeled. The duty cycle forM. ^ 10'^ Mq 
is - 0.2 at z 4 - 5, falling to 0.03-0.08 at z = 2 - 3, when 
quasar activity is at its peak, then dropping to 0.003 at z = 1 
and ^ lO^** at z = 0. Below z 3, the duty cycle rises to- 
wards low black hole masses. This "downsizing" evolution, 
in which high mass black holes complete their growth early 
but low mass black holes continue to grow at late times, is 
required by the observed LF evolution in any model with 
approximately constant mo. Above lO^M© the duty cycle 
curves are generally flat, because the constant slope of the 
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bright end of the LF drives the growth of a mass function 
with a parallel slope. Below z = 1, the LF slope changes, 
and the duty cycle curves are no longer flat. By running sev- 
eral cases with varying initial conditions, we find that the duty 
cycles for z < 3.5 are insensitive to the assumed initial duty 
cycle (which determines the initial mass function) provided 
0.01 <Po(z = 6) < 1. 

Several lines of independent observational evidence have 
set constraints on the duty cycle of AGNs at different red- 
shifts. Porciani et al. (2004), Croom et al. (2005), Porciani 
and Norberg (2006), and da Angela et al. (2006) have ana- 
lyzed large samples of optical AGNs in the Two Degree Field 
(2dF) QSO Redshift Survey, showing that their large-scale 
(> 1 Mpc) clustering properties are consistent with quasars 
residing in very massive halos withMj, ^5x10'^— lO'^M©. 
Comparing the AGN abundance to the abundance of such ha- 
los implies duty cycles of order of ^ (1 —4)% in the red- 
shift range z ^ 1 — 2. Similar results have been confirmed by 
large-scale optical AGN clustering studies in the Sloan Dig- 
ital Sky Survey (SDSS; e.g., Myers et al. 2007). Shen et al. 
(2007) have extended the SDSS results to z = 3 - 4, obtain- 
ing similar results for halo masses and duty cycles. Similar 
or lower values for the duty cycle could be found directly by 
comparing the AGN number density to the galaxy luminos- 
ity function at comparable redshifts, the latter now well es- 
timated and complete down to rather low luminosities (e.g., 
Pannella et al. 2006). Using Chandra X-ray observations of 
the Hubble Field North, Nandra et al. (2002) and Steidel et 
al. (2002) find that about 3% of a large spectroscopic survey 
of ^ 1000 Lyman Break Galaxies at the average redshift of 
z ~ 3 host an X-ray bright, moderately obscured AGN. Within 
the uncertainties of the comparison, these estimates of AGN 
duty cycles are in good agreement with the predictions in Fig- 
urelTJ;. We will present predictions for AGN clustering bias in 
§ ID below, and we will carry out quantitative comparisons to 
observed clustering in future work. The occurrence of AGN 
in massive star-forming galaxies is much higher than these 
overall duty cycle estimates (Alexander et al. 2003; Borys et 
al. 2005), supporting the scenario in which star-formation in 
massive galaxies is closely related to black hole growth (e.g., 
Granato et al. 2006). 

Figure|7}l shows the black hole accretion histories as a func- 
tion of relic mass and redshift, following equation (fT2l) . Trac- 
ing back one of these curves shows the mass that a "typical" 
black hole of present mass M, had at earlier times, as pre- 
dicted by our model. The dot-dashed line shows the minimum 
black hole mass associated with Lmin(z) at each redshift. By 
construction, all growth curves are flat (dashed lines) before 
crossing Lmin(z)- Figure|7}i shows again that high mass black 
holes build their mass early and experience little growth at 
late times, while lower mass black holes grow rapidly at lower 
redshifts. 

4.3. Parameter variations 

The top panels of Figure [8] show the effect of varying the 
two main model parameters, the radiative efficiency e and ac- 
cretion rate riid. Figure [8^ compares models with varying e 
and accretion rate fixed to the reference model value rh — 0.60. 
Lowering (raising) the radiative efficiency has the effect of 
shifting the accreted mass function up (down) and shifting the 
peak of the mass distribution to higher (lower) masses. Re- 
call that the luminosity Eddington ratio is A = m/o,i; if we 
held A fixed instead of rii, then curves in Figure [8^ would be 
vertically displaced with no horizontal shift. Figure[8}5 shows 



models with varying mo, all with e = 0.065. The Soltan ar- 
gument impHes that p, = J M,^{M,)dlogM, should be the 
same for models with the same e that reproduce the same ob- 
served LF. The area under the curves in Figure [Sj) is thus the 
same for all models, but the peak in the mass function shifts 
to lower M, for higher otq because a given observed luminos- 
ity is associated with a lower black hole mass. While there are 
still systematic uncertainties in the normalization and shape of 
the local mass function, the peak of the M,^{M, ) distribution 
(Figure |5]) is in good agreement among all the estimates, and 
this, in turn, sets strong constraints on mo. Shifting the peak 
of the accreted mass distribution to log(M,/MQ) ^ 8 would 
require Super-Eddington accretion rates, while shifting it to 
log(M,/MQ) ^ 9 requires mo ~ 0.2 — 0.3. Changes to e and 
niQ shift the accreted mass function in different directions in 
the M,<i>(M, ) — M, plane, but they do not alter its shape. 

The lower panels of Figure [8] show the effect of chang- 
ing the input LF. In panel [8]:, solid and dotted curves show 
models with no Compton-thick sources or a double fraction 
of Compton-thick sources, respectively, with radiative effi- 
ciency and accretion rate fixed at the reference model values 
e = 0.065 and mo = 0.6. The predicted black hole mass func- 
tions are similar in shape but offset in amplitude by about ± 
20%. These offsets are similar to our estimated uncertainty 
in the local mass function, shown by the grey band. We can 
restore agreement between the double Compton-thick model 
and the central estimate of the local mass function by raising e 
to 0.08 (to produce more luminosity with the same black hole 
mass) and lowering nio to 0.5 (to compensate the influence of 
higher e on the peak location). Conversely, the no-Compton- 
thick model yields similar predictions if e is lowered to 0.05 
and mo is raised to 0.8. As noted in §|2l we find better agree- 
ment with the observed X-ray background for no Compton- 
thick and worse agreement for double Compton-thick (see 
Figure O. 

Figure [8}J compares our reference model to the one ob- 
tained by integrating the AGN bolometric LF of HRH07. For 
the reference model values e = 0.065 and mo ~ 0.6, the pre- 
dicted mass function shifts to higher normalization and higher 
peak mass (solid line). Principally this difference reflects the 
higher bolometric correction adopted by HRH07, which shifts 
their LF to higher normalization at 1 < z < 3 (Figure HI . The 
HRH07 estimate also has higher peak luminosity in this red- 
shift range, leading to higher peak mass. Adopting a higher 
efficiency and accretion rate, e = 0.09 and mo = 1, largely 
compensates these two effects. As noted in § |2l the HRH07 
LF has a steeper bright-end slope at intermediate redshifts, 
and this leads to a steeper high-mass slope of the accreted 
mass function, improving agreement with observational esti- 
mates above M, ^ lO^M©. 

To present our model comparison in a more global way, we 
characterize the local mass function by four quantities that 
contain complementary information: 

• The integrated black hole mass density p, = 
/ M,^{M,)dlogM,. 

• The location logMpEAK at which the mass function 
peaks in the M,$(M,) -M, plane. 

• The width of the mass function, characterized 
by ( A log M,) such that twice the integral 
from logMpEAK to logMpEAK + (AlogM,),/2 
contains half the total mass density, 2 x 
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/logMpL HM,)M,d\ogM,= \/2 X p,. If M.<i>(M.) 

were a Gaussian, (AlogM,)i/2 would be ~ 0.26 times 
the full width at half maximum. 

• The asymmetry of the mass function, characterized by 
Ap,/p,, the normalized difference in the mass density 
integrated above and below logMpEAK- 

Figure |9] plots the predicted values of these four quantities 
as a function of the radiative efficiency e for models with our 
standard AGN LF and accretion rates nio = 1,0.6,0.45,0.3. 
Grey horizontal bands show the observational estimates corre- 
sponding to the grey band in Figure|5] Short-dashed and long- 
dashed horizontal lines show instead the parameters for the lo- 
cal mass functions inferred by Graham et al. (2007) and Hop- 
kins et al. (2007b), respectively (see Figure |5]l. As expected 
from the Soltan argument, the predicted p, is proportional to 
(1 — e)/e but independent of otq- Consistency with the grey 
band observational estimate requires 0.063 < e < 0.1.^ Re- 
producing the observed value of logMpEAK/M© k, 8.5 then 
requires 0.4 < mo < 0.8, though for any given e the value of 
m^) is determined to ^ 10%. As noted in our discussion of Fig- 
ure [8^, the location of the mass function peak is determined 
by A = L/LEdd largely (but not completely) independent of e, 
and reproducing the observed logMpEAK implies A « 0.4 — 0.5 
over the range 0.05 ^ e < 0.13. 

The higher p, implied by the Hopkins et al. (2007b) lo- 
cal mass function requires a lower efficiency, e w 0.06 which 
in turn implies a higher my w 1 to match logMpEAK. (Note 
that this is not the comparison in Figure |8}l, where we use the 
HRH07 luminosity function but show our standard estimate of 
the local mass function). Our inferred ranges for e and mo are 
consistent with the findings of Marconi et al. (2004) and S04, 
who obtained (e, A)=(0.08, 0.5) and (0.09, 0.3), respectively. 
Note that S04 tends towards slightly higher values of the ra- 
diative efficiency mainly because they used an X-ray bolomet- 
ric correction normalized to the optical bolometric correction 
of Elvis et al. (1994), i.e., Cb=11.8, while we have adopted 
the lower value suggested by recent work (see §|2]l. 

With e and mo fixed by matching p, and logMpEAK, there 
is no freedom in the model to adjust the predicted width and 
asymmetry of the mass function, but these prove remarkably 
insensitive to e and nio in any case. To substantially alter these 
predictions one must either change the input luminosity func- 
tion or add a new physical ingredient to the model, such as 
mass- or redshift-dependent mo, mergers, or multiple accre- 
tion rates. For our minimal model and standard LF estimate, 
the predicted widths are larger than any of the observational 
estimates. This discrepancy is largely a consequence of the 
shallow high-mass slope of the predicted mass function (see 
FigurelH), which in turn is sensitive to the bright-end slope of 
the luminosity function in the range 1 < z < 3, as discussed 
earlier. Note, however, that the observational estimates of 
(AlogM,)i/2 depend on the extrapolation of the black hole- 
galaxy correlations above M, « 10^ M0, where they are quite 
uncertain, so the grey band in Figure |9] should not necessarily 
be taken as a true upper limit. The predicted asymmetries are 
in reasonable agreement with the observational estimates, ex- 
cept for the Graham et al. (2007) mass function, in which the 
population of low mass black holes is greatly reduced. 

' Our reference model of j) l4.2l has e at the low end of this range, rather 
than the middle, because we chose parameters to match the mass function 
near its peak, where it is best determined. However, the reference model 
predicts p. = 5.24 X 10^ Mq Mpc~^, above the middle of the grey band. 



Figure [TO] presents a similar comparison for models with 
different input luminosity functions. Here we fix mo — 0.75 
in all cases, noting that the value of mo affects the predicted 
logMpEAK but has little impact on other quantities. Elimi- 
nating Compton-thick sources from our standard LF (dotted 
curve) lowers the efficiency required to match a given p, 
by about 20%, to e = 0.062 for p, = 4.26 x lO^MoMpc-^ 
(short-dashed line). Doubling the Compton-thick contribu- 
tion (not shown) has a similar effect in the opposite direction, 
raising the required e by 20%. Changing the Compton-thick 
contribution has little impact on logMpEAK for a given e, or 
on the width or asymmetry of the predicted mass function. 
Eliminating all obscured sources and retaining only Type I 
quasars drastically reduces the predicted p,, so such a model 
(which is implausible on direct observational grounds any- 
way) can only be reconciled with the local black hole popu- 
lation if the efficiency is very low or if the true value of p, 
is significantly below our observational estimate. However, a 
model in which the full bolometric luminosity function is just 
three times that of Type I AGN — two obscured sources for 
every unobscured source independent of luminosity and red- 
shift — yields similar predictions to our full model with the 
luminosity-dependent column density distribution of U03. 

For the HRH07 luminosity function, we require 0.079 < 
e < 0. 125 to reproduce the grey-band range of p., 25% higher 
than our standard LF at fixed p.. As discussed earlier, this 
difference mainly reflects the higher bolometric correction 
adopted by HRH07, which boosts the normalization of the 
LF at intermediate redshifts. The HRH07 LF also leads to a 
higher MpEAK at a given e, and matching logMpEAK/A^o ~ 8.5 
implies 0.7 ^mo 1.1 for the range of e that matches p.. 
The HRH07 LF predicts a somewhat narrower and more 
asymmetric mass function, mainly because of the steeper 
bright-end slope at intermediate redshifts, which steepens the 
high mass end of the black hole mass function. These re- 
sults accord with those shown already in Figure [8}J. If we 
adopt both the HRH07 luminosity function and the Hop- 
kins et al. (2007b) black hole mass function, which has 
p, = 5.7 X lO^MoMpc"^ and logMpEAK/M© « 8.4, we find 
e = 0.075 and mo = 1.3. 

The remaining uncertainties in the luminosity function and 
the local black hole mass function still leave considerable un- 
certainty in the radiative efficiency. If we assume that our 
no-Compton-thick LF and the HRH07 LF bracket the lumi- 
nosity function uncertainty and our grey band brackets the 
p, uncertainty, the allowed range is 0.05 < e < 0.125, and 
there is clearly some room for the luminosity function or p, 
to be outside this range. However, the high efficiencies pre- 
dicted by MHD simulations of thick accretion disks around 
rapidly spinning black holes, e « 0.15 — 0.2 (e.g., Gammie 
et al. 2004), can only represent the mean efficiency of black 
hole accretion if p, is considerably below our estimates or 
the luminosity function is substantially higher than even the 
HRH07 estimate. There are some p, estimates as low as 
2 - 3 X IO^MqMpc-3 (e.g., Wyithe et al. 2004). Our con- 
clusions about efficiency are in agreement with most previous 
studies, though Elvis et al. (2002) reached a different conclu- 
sion mainly because they assumed that the sources produc- 
ing the X-ray background have an effective mean redshift of 
z^2, while subsequent data have shown that the peak contri- 
bution to the background is at z < 1 (e.g., U03; La Franca et 
al. 2005). 

4.4. Redshift-dependent Eddington ratio 
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So far we have assumed that all black holes at all times ra- 
diate at a constant, mildly sub-Eddington rate, an assumption 
supported by Kollmeier et al. (2005)'s results for luminous 
quasars in the AGES survey at 0.5 ^ z ^ 4. However, Edding- 
ton ratios of local AGN have a wide distribution, with a large 
fraction of Seyferts radiating at A < 0.1. Some studies also 
show evidence for mass-dependence of the typical Eddington 
ratio (e.g., Heckmann et al. 2004; Vestergaard 2004; McLure 
& Dunlop 2004; Constantin & Vogeley 2006; Dasyra et al. 
2007; Shen et al. 2007b), though systematic uncertainties in 
the reverberation mapping techniques and extrapolations of 
empirical virial relations (e.g., Kaspi et al. 2000; Bentz et 
al. 2006) make it hard to draw firm conclusions. In this sec- 
tion and the one that follows, we discuss simple models with 
mo values (and thus Eddington ratios) that depend on redshift 
or mass. More realistic models would incorporate a range 
of Eddington ratios with the relative importance of high and 
low accretion rates depending on mass or redshift; we will in- 
vestigate such models in future work. We note again that the 
integrated black hole mass density should depend only on e 
(equation HI, but the shape of the black hole mass function 
can change if nii) is not constant. 

Figure [TT| compares the results of our reference model, 
which has e ~ 0.065, niQ ~ 0.6, to a model with the same e 
but 



m{z) = [1 -exp(-z/zs)], 



(17) 



with Zs = 2. This implies values of niQ = 

0.78,0.63,0.39,0.22,0.095,0.049 at z = 3,2,1,0.5,0.2,0.1, 
respectively, which are consistent with the average drop 
shown in Figure 6b of Vestergaard et al. (2004) for the 
more massive black holes, assuming the distribution starts 
at ;«o ~ 1 at very high redshifts. The decreasing mo{z) 
significantly alters the shape of the black hole mass function 
(Figure [TTh) by shifting the low redshift growth of the mass 
density away from low mass black holes and towards higher 
mass black holes. This model still requires a substantial 
number of low-mass black holes residing in spiral galaxy 
bulges if it is to match the local mass function. 

Figures [TTb compares the duty cycles produced by the 
in(){z) model to those from the reference model, at several 
redshifts. The reference model duty cycles are strongly de- 
creasing functions of redshift and black hole mass, dropping 
by up to two orders of magnitude in the lower redshift bins. 
The mo (z) -model "damps" the downsizing, leaving a much 
milder variation of the duty cycle over time and mass. As 
emphasized by SW03, the downward shift of the LF break 
luminosity at z < 2 can be explained either by a strong sup- 
pression of activity in high mass black holes or by a drop in 
characteristic Eddington ratios. Our mo (z) model is a compro- 
mise between these two explanations; higher mass black holes 
still build their mass earlier than low mass black holes, but the 
difference is smaller than in the reference model (FigurefTTb). 

Greene & Ho (2007) have recently estimated the local MF 
for active black holes through the virial relations mentioned 
above. When compared with the local mass function of relic 
black holes, their study supports a duty-cycle of Pq ~ 10~^^^ 
for black holes with mass logM,/M0 > 7, consistent with 
the z ^ duty-cycle inferred from this model, and about two 
orders of magnitude above the one corresponding to our ref- 
erence model (FigurefTTb). 

4.5. Eddington ratio varying with black hole mass 



In the previous section we studied a model in which the 
Eddington ratio steadily decreases with redshift and is fixed 
with black hole mass. Here we consider the complementary 
model in which the Eddington ratio distribution is constant 
with redshift but declines with black hole mass. 

We adopt 

mo(M.) = 0.445 x(^y' (18) 

with a constant radiative efficiency e = 0.075. As discussed 
in § 14. 1[ we evolve this model by integrating equation ( fTSl l 
with a fixed black hole mass grid, extrapolating the luminosity 
function below Lmin(z). Relative to the reference model, the 
m^){M,) model associates the peak of the LF with higher mass 
black holes, especially at lower redshifts. The low redshift 
growth of low mass black holes is strongly suppressed, and 
the z = abundance of low mass black holes is much lower, 
consistent with the Graham et al. (2007) local mass function 
estimate, as shown in Figure [TJ] The dependence of duty 
cycle on mass is also much stronger in this model because of 
the paucity of low mass black holes; note that the duty cycle at 
low masses is high but the amount of growth is low because 
of the low accretion rates. Getting our mass-dependent mo 
model to run at all requires a high initial black hole density, 
with a duty cycle Po(z = 6) = 0.01 — 0.02; otherwise, the slow 
growth of low mass black holes leads to required duty cycles 
that exceed unity at intermediate redshifts. Alternatively, we 
can start with a higher initial duty cycle and only "turn on" the 
mass dependence of equation ( fTSt after an epoch of growth at 
constant niQ. 

More generally, we can vary the Eddington ratio both in 
redshift and mass, mo(M,,z) = T{z)B{M,). The continuity 
equation can then be written 

dn{M,,t) _ 



d 



dlogM, 
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where we have set to zero the second derivative 
assuming a power-law relation like 
that in equation (fTsT l. We have used equation (19) to inves- 
tigate models in which we vary the Eddington ratio with 
redshift and mass as given in equations ( [TtI i and ( fTSl l. We 
find that we can get somewhat better fits to Graham et al.'s 
(2007) local mass function, since both the implemented 
trends act to decrease the numbers of low mass black holes. 

4.6. Black hole mergers 

Observations show, and hierarchical galaxy formation mod- 
els predict, that a significant fraction of galaxies experience 
mergers with comparably massive galaxies during their life- 
time. At least some of these galaxy mergers are likely to be 
accompanied by mergers of the central black holes that they 
contain, though the mechanisms that shrink black hole or- 
bits to the scale where gravitational radiation can drive a final 
merger are not fully understood (see Merritt & Milosavlje- 
vic 2005). In a future paper, we will will combine our ac- 
cretion model with theoretically predicted merger rates for 
cold dark matter subhalos (Yoo et al. 2007) to create mod- 
els with realistic contributions of accretion and merger driven 
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growth. Here we illustrate the potential impact of mergers 
on the black hole mass function using a simple mathematical 
model that assumes constant probability of equal mass merg- 
ers per Hubble time, similar to the models of Richstone et 
al. (1998) and SW03. We will ignore physically interesting 
complications such as ejection of black holes by gravitational 
radiation or three-body interactions or variations of radiative 
efficiency caused by the impact of mergers on the black hole 
spin distribution (e.g., Hughes & Blandford 2003; Volonteri 
et al. 2003; Gammie et al. 2004; Islam et al. 2004; Yoo & 
Miralda-Escude 2004; Merritt & Milosavljevic 2005; Volon- 
teri et al. 2006). These effects, along with unequal mergers 
and the mass and redshift-dependence of merger rates, should 
be considered in a complete model of the evolving black hole 
population. Mergers redistribute mass within the black hole 
population, but they do not change the integrated mass den- 
sity, so they do not affect the integrated Soltan (1982) argu- 
ment. (In principle, gravitational radiation during mergers can 
reduce the integrated mass density [see Yu & Tremaine 2002], 
but we do not consider this effect here). 

In our model, a black hole of mass M, has a probability 
Anerg of merging with an equal mass black hole in the Hubble 
time tniz) (age of the Universe at redshift z). Therefore the 
fraction F of black holes that merge in a timestep A? is 

Af 

F = i^merg X — - . (20) 

At each time fi, we first advance the mass function to time t2 = 
f 1 + Af with accretion only, then add to each bin of the mass 
function an increment (which may be positive or negative) 

A$(M.,f2) = ^ ' ' ^ -Fx<^{M„t2), (21) 

where the second term represents black holes lost from the bin 
by merging and <1> (^y- , f2) is calculated by interpolation. 

Figure [13] shows the evolution of the black hole mass func- 
tion for a model with Pmeig — 0.5 and our reference model val- 
ues of e = 0.065 and mo = 0.6. Comparing the z — 0.02 out- 
put to that of the reference model shows that the net effect of 
merging is to slightly lower the abundance of black holes be- 
low the peak of the mass function and to significantly increase 
the abundance of very massive black holes, as expected (see 
also Malbon et al. 2006; Yoo et al. 2007). Although we start 
merging at z = 6, comparison of Figure[T3]to FigurelTJ) shows 
that the mass function is nearly identical to that of the refer- 
ence model down to z = 1 . Mergers at this level only change 
the mass function noticeably once accretion-driven evolution 
has slowed. 

Since our reference model already produces an excess of 
massive black holes relative to local estimates of the black 
hole function, adding mergers only makes the match to ob- 
servations worse. However, the impact of mergers is evident 
mainly at M, > 10''Mq, where the local mass function esti- 
mates are most sensitive to the adopted scatter in the black 
hole-host scaling relations and to the extrapolations of these 
relations to the most luminous galaxies. Accommodating the 
predictions of our merger model would require an intrinsic 
scatter of ^ 0.5 dex at high masses, or a change in slope of 
the scaling relations. Either of these changes seems possi- 
ble given the current observational uncertainties, but neither 
seems likely (see, e.g., discussions by Wyithe 2004; Batchel- 
dor et al. 2006; Lauer et al. 2006; Tundo et al. 2007). If we 
adopt the HRH07 LF in place of our standard LF, then the pre- 
dicted black hole mass function with mergers is similar to that 



of our reference model without mergers (see Figure [T3]). The 
impact of mergers at this level is therefore comparable to the 
systematic uncertainties associated with the AGN luminosity 
function. 

In fact, the value Pmerg = 0.5 is high compared to theo- 
retical predictions for massive galaxies (Mailer et al. 2006), 
and theoretically predicted merger rates decline towards lower 
masses. Observational estimates of the galaxy merger rate and 
its mass dependence span a substantial range (see, e.g.. Bell 
et al. 2006; Conroy et al. 2007; Masjedi et al. 2008, and 
references therein), and Pmerg = 0.5 is roughly consistent with 
the high end of these estimates. We conclude that the im- 
pact of mergers on the black hole mass function is probably 
small compared to remaining uncertainties in accretion-driven 
growth, except perhaps for the rare, high mass black holes. 
The most interesting impact of black hole mergers may arise 
indirectly, through their effect on black hole spins and thus on 
radiative efficiencies (e.g., Volonteri et al. 2005). 

4.7. Tabulation of luminosity and mass functions 

For the convenience of readers who may wish to use them, 
we provide electronic tables that list our estimate of the 
AGN bolometric LF and some of our model predictions of 
the black hole mass function and duty cycle, all as a func- 
tion of redshift. Table |2] lists our standard LF estimate, the 
LF after eliminating or doubling the number of Compton 
thick sources, and the HRH07 LF estimate evaluated at the 
same redshifts and luminosities for convenient comparison. 
Table [3] lists the predicted black hole mass function at the 
same redshifts for our reference model values of e = 0.065, 
riiQ ~ 0.60, computed for each of the luminosity functions in 
Table |2] The final column lists the mass function predicted 
for the HRH07 LF and the parameter values e — 0.09 and 
riiQ — \, which yield a good match to the local mass func- 
tion for this LF. Table |4] lists instead the duty cycles corre- 
sponding to the same models reported in Table [5] computed 
at the same redshifts and black hole masses. [Prior to pub- 
lication, the full tables can be found in electronic format at 
http : / / www . astronomy . ohio-state . edu/ 
'^shankar /Models.] 

5. SPACE DENSITY AND BIAS OF AGN HOSTS 

Predictions of the black hole mass function at redshifts 
z > can be tested observationally if one assumes that black 
hole mass increases monotonically with the luminosity of the 
host galaxy or the mass of its parent halo. This assumption 
is unlikely to be perfect, but given the tight correlation be- 
tween black hole mass and bulge mass observed today, it may 
be a reasonable approximation. Figure [l4h shows, at several 
redshifts, the predictions of our reference model for the cumu- 
lative space density of galaxies that can host an AGN of lumi- 
nosity L or greater Since the model assumes a single mo, this 
is simply equal to the cumulative space density of black holes 
of massM. >L/(/o.i;7ioO = 1-9 x 10** Mq x (L/lO^^ergs-i). 
If the monotonic assumption holds, then the observed space 
density of galaxies brighter than Lhost(i'AGN) should equal the 
space density predicted in Figure [l4h. where Lhost(i'AGN) is 
the luminosity of galaxies that host AGN of luminosity Lagn- 

As emphasized by Haiman & Hui (2001) and Martini & 
Weinberg (2001), the clustering of AGN can be a powerful 
diagnostic for the duty cycle of black hole activity: a high 
duty cycle implies that halos of AGN are rare, hence high 
mass, hence strongly clustered. In the context of our models, 
we can predict the halo mass M/, associated with black holes 
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of mass M, by matching the cumulative mass functions, 

$(>M.,z)-$(>M,„z). (22) 

Here $(> M,,z) is the space density of black holes more mas- 
sive than M, at redshift z, in units of Mpc^^, as predicted by 
our evolutionary model, and $(> Mk,z) is the space density 
of halos more massive than Mh expected for a ACDM cos- 
mological model, which we compute using the Sheth & Tor- 
men (1999) halo mass function with cosmological parameters 
Vl„, — 0.3, f^A = 0.7,/; = 0.7, (Ts = 0.8 and linear power spec- 
ti-um taken from Smith et al. (2003) with T = 0.2. This de- 
termination assumes both a monotonic relation between black 
hole mass and halo mass and one black hole per halo, but we 
have checked that allowing black holes to reside in subhalos 
does not substantially change our results if we adopt the sub- 
halo statistics of Vale & Ostriker (2004). Finally, we compute 
the bias of black holes of mass M,, and thus the bias of AGN 
of luminosity L — fo imolM,, using Sheth et al.'s (2001) ana- 
lytic model for halo bias, 

b{L,z)=b{M.,z)=b{Mh,z). (23) 

Figure fT4b shows the predicted bias as a function of luminos- 
ity L at several redshifts for the reference model of § 14.21 The 
predicted bias is much stronger at high redshift because the 
comoving space density of black holes is lower (FigurelT) and 
because the bias of halos of a given space density is higher 
Points in Figure [T4b show recent observational estimates of 
AGN clustering bias from the 2dF and SDSS quasar redshift 
surveys (da Angela et al. 2007; Myers et al. 2007; Por- 
ciani & Norberg 2007) in the redshift range 0.8 ^z^l3 and 
1-7 < z < 2.1 (Porciani & Norberg 2007) and from SDSS at 
3 < z < 3.5 (Shen et al. 2007a). Our reference model is quite 
successful at predicting the absolute level of bias (and thus 
AGN clustering strength) for quasars with L > 10'*^ ergs"' 
at z ^ 1, z ^ 2 and z ^ 3, and, especially, at explaining the 
strong trend of bias among these samples. It appears to un- 
derpredict the clustering of lower luminosity AGN at z ^ 1, 
though Myers et al. (2007) note the including the possible 
impact of stellar contamination on their sample expands their 
error bars by a factor of 1 .5 (an effect we have not included 
in Figure [l4b). 

We reserve a detailed examination of AGN clustering con- 
straints — including an assessment of what models can be 
ruled out with present data or improved measurements — for 
future work. Allowing scatter between black hole mass and 
halo mass or a range of Eddington ratios can significantly alter 
predictions both for AGN host density and for AGN cluster- 
ing, so these classes of models will be especially interesting 
to explore. For example, the stronger clustering in the Myers 
et al. (2007) sample relative to our reference model predic- 
tions could indicate that a fraction of these lower luminosity 
AGN are associated with high mass black holes radiating at 
low Eddington ratios (e.g., Lidz et al. 2006). 

6. DISCUSSION 

We have constructed self-consistent models for the evolu- 
tion of the supermassive black hole population and the AGN 
population, in which black holes grow at the rate implied 
by the observed luminosity function given assumed values 
of the radiative efficiency e and the characteristic accretion 
rate mo = M,/MEdd (see Eqs. [TO] [121 [TsT l. These models 
can be tested against the mass distribution of black holes in 
the local universe, and they make predictions for the duty cy- 
cles of black holes as a function of redshift and mass, which 



can be tested against observations of quasar hosts and quasar 
clustering if one assumes an approximately monotonic rela- 
tion between the masses of black holes and the masses of 
their host galaxies and halos. Our method is similar to that 
used previously by Cavaliere et al. (1971), Small & Bland- 
ford (1992), Yu & Tremaine (2002), SW03, Marconi et al. 
(2004), and S04. However, we have drawn on more recent 
data on the AGN luminosity function and the local black hole 
mass function, and we have considered a broader spectrum 
of models. This approach can be considered a "differential" 
generalization of the Soltan (1982) argument relating the in- 
tegrated emissivity of the quasar population to the integrated 
mass density of the local black hole population. 

Our model for the bolometric AGN luminosity function 
starts from the model of U03, based on X-ray data from a 
variety of surveys, but we adjust its parameters as a func- 
tion of redshift in light of more recent measurements and 
data at other wavelengths. In agreement with previous stud- 
ies we find that the LF of optical AGNs is roughly consistent 
with that of X-ray AGN that have absorbing column densities 
logA^///cm^^ < 22 and that unobscured AGN dominate the 
bright end of the LF. We show that the latest constraints on the 
hard X-ray background (£ - 10 - 100 keV) from INTEGRAL 
and from the PDS instrument on BeppoSax support a reduced 
normalization relative to extrapolations from other missions 
at lower energies. They therefore favor a lower contribution 
from very highly obscured AGN (logA^///cm^^ > 24.5) than 
some previous estimates (but see also Gilli et al. 2007). Our 
estimate of the bolometric AGN LF is independent in imple- 
mentation but similar in spirit to that of HRH07. The most 
significant differences for purposes of this investigation are 
that HRH07 have a higher LF normalization at z 1 — 3, prin- 
cipally because of their choice of bolometric correction, and 
they have a steeper bright-end slope at z 1—2.5, where they 
adopt the slopes measured by Richards et al. (2006) from the 
SDSS and we use the slopes inferred by U03 from X-ray data. 
We regard the differences between our estimate and that of 
HRH07 as a reasonable indication of the remaining system- 
atic uncertainties in the bolometric LF of AGNs. 

With our LF estimate, the bolometric emissivity of the AGN 
population tracks recent estimates of the cosmic star forma- 
tion rate as a function of redshift. (Comparisons based on 
the space density of high luminosity quasars [e.g., Richards et 
al. 2006, Osmer 2004 and references therein] reach a different 
conclusion because at low redshifts the bright end of the AGN 
LF drops much more rapidly with time than the overall emis- 
sivity.) The integrated black hole mass density implied by this 
emissivity is ^ 8 x lO^'* of the stellar mass at all redshifts, or 
about half of the estimated ratio of black hole mass to bulge 
stellar mass in local galaxies. This tracking favors scenarios 
in which black holes and the stellar mass of bulges grow in 
parallel, with about 50% of the star formation linked to black 
hole growth at all redshifts. These findings are hard to rec- 
oncile with any models where black hole growth substantially 
precedes stellar mass buildup, or with recent claims that the 
ratio of black hole mass to stellar mass is much larger at high 
redshift than the local value. However, our finding refers to 
integrated densities, so it does not indicate the relative timing 
of black hole and bulge growth on an object-by-object basis. 

Observational estimates of the local black hole mass func- 
tion still show substantial discrepancies among different au- 
thors, depending on the correlation used to derive it (e.g., 
M, — a, M, — Lbuige, M, — Mstar, M,-Sersic index, or funda- 
mental plane), the calibration of the correlation, and the intrin- 
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sic scatter of the correlation. Above M, ^ 1O''M0, estimates 
depend on extrapolation of the observed scaling relations, and 
below M, ~ IO^'^Mq they are sensitive to the treatment of 
spiral bulges. The grey band in Figure|5](and subsequent fig- 
ures) encompasses most estimates, but the fundamental plane 
(Hopkins et al. 2007b) and Sersic index (Graham et al. 2007) 
methods imply more sharply peaked mass functions. The inte- 
grated mass densities of all of these estimates are in the range 
p. ~ 3 -5.5 X lO^MoMpc-l 

Our simplest models assume a single characteristic Edding- 
ton accretion rate mo, independent of redshift and black hole 
mass, and all of our models assume a single value of the 
radiative efficiency e. Matching the local black hole mass 
density requires e = 0.075 x {p,/4.5 x IO^MqMpc"^)"' for 
our standard estimate of the AGN LF, or e ~ 0.094(/9,/4.5 x 
10-**MqMpc-3)-i for the HRH07 luminosity function.^ With 
e thus fixed, the value of mo determines the peak of the pre- 
dicted local black hole mass function in the M,<I>(M,) vs M, 
plane. Note that with our definitions the Eddington luminosity 
ratio is A = L/L-Edd ~ ;7Jo(e/0.1) (equation |9]l. Matching the 
observed peak at logM,/MQ ^ 8.5 implies rito « 0.6 (A « 
0.45) for our standard LF estimate and mo « 1 (A ~ 0.95) for 
the HRH07 LF. Lower values, mo = 0.1 — 0.3, shift the peak 
location to untenably high masses, logM,/M0 ^ 8.9 — 9.3 or, 
for HRH07, logM./Mo - 9.1-9.6. The single-mo models 
achieve a reasonable match to our "grey-band" observational 
estimates of the width and asymmetry of the local mass func- 
tion, though for our LF estimate the predicted mass function is 
too high at M, > 10'' M0 and is therefore somewhat too broad. 
Single-;7io models cannot reproduce the more sharply peaked 
local mass functions estimated by Graham et al. (2007) or 
Hopkins et al. (2007b). 

Our reference model, which has e = 0.065, mo = 0.60, and 
our standard LF estimate, predicts a duty cycle for activity 
of 1O''M0 black holes that declines steadily from 0.15 at 
z = 4 to 0.07 (z = 3), 0.035 (z = 2), 0.004 (z = 1), and lO^^ 
(z = 0). The decline in duty cycle for lower mass black holes 
is much shallower. Massive black holes therefore build their 
mass relatively early while low mass black holes grow later, 
the phenomenon often referred to as "downsizing". Our re- 
sults on mean radiative efficiency and duty cycle evolution 
are also in qualitative agreement with those found by Haiman 
et al. (2004). The predicted duty cycles seem in reasonable 
accord with observational estimates, though these estimates 
have considerable uncertainty and do not, as yet, probe mass 
and redshift dependence in much detail. The electronic tables 
described in § 14.71 provide tabulations as a function of red- 
shift of our AGN bolometric LF estimate, the HRH07 LF, and 
black hole mass functions and duty cycles for single-jiio mod- 
els that are in good agreement with the observed z = mass 
functions given these LF inputs. 

We have examined models in which the Eddington ratio ac- 
cretion rate mo is reduced at low redshift or at low black hole 
mass. Declining redshift evolution of mo damps "downsiz- 
ing," reducing the dependence of duty cycles on black hole 
mass and redshift. This model produces a typical duty cycle 
Pq ~ 10^^'^ at z = 0, about two orders of magnitude higher 
than in the reference model and consistent with the local duty 
cycle inferred from observations by Greene & Ho (2007). 
In general terms, the observed luminosity-dependent density 
evolution of the AGN LF can be explained by preferential sup- 

* More precisely, it is — e) that is proportional to p7' > but tlie differ- 
ence from e oc p^' is tiny over tlie allowed range. 



pression of activity in high mass black holes at low redshift, 
by a decline in the typical accretion rate at low redshift, or 
by some combination thereof. The mass-dependent riio model 
associates more of the AGN emissivity to high mass black 
holes, so it predicts a z = mass function that is more sharply 
peaked, in better agreement with the estimates of Hopkins et 
al. (2007b) and Graham et al. (2007) but worse agreement 
with other estimates. This model predicts a stronger mass de- 
pendence of duty cycles than our reference model because it 
maps low mass black holes, whose abundance is already sup- 
pressed, to less luminous, and hence more common, AGN. 

We have also considered a model in which each black hole 
has a 50% probability per Hubble time of merging with an- 
other black hole of equal mass. Merger-driven growth in this 
model has little impact on the black hole mass function until 
z < 1, when accretion-driven evolution has slowed. Low red- 
shift mergers slightly depress the low mass end of the z = 
mass function and significantly enhance the high mass tail, 
worsening the agreement with observational estimates. Mod- 
els incorporating theoretically predicted merger rates can al- 
low more realistic calculations of the impact of mergers on the 
black hole population; this impact will probably be smaller 
than in the simplified model considered here. 

We have calculated the clustering bias of AGN as a function 
of luminosity and redshift for our reference model, assuming 
a monotonic relation between black hole mass and halo mass. 
The predictions are in reasonable accord with observational 
estimates. We will examine AGN clustering predictions in 
more detail in future work, with attention to what models can 
be excluded by the data and what can be learned by matching 
the full AGN correlation function in addition to an overall bias 
factor 

MHD simulations (e.g., Gammie et al. 2004; Shapiro 
2005) show that disk accretion onto Kerr black holes spins 
them up to an equilibrium spin rate a w 0.95 (where a ~ \ 
is the angular momentum parameter for a maximally rotat- 
ing black hole). The radiative efficiency in these models is 
e « 0.16 — 0.2. These high efficiencies would lead to black 
hole mass densities a factor of two or more below our central 
estimate, and below our estimated lower bound. Furthermore, 
our results show that models with e in the range 0.06 — 0.1 1 
can achieve a good match to the overall shape of the black 
hole mass function near the peak in M,$(M,), not just the 
value of p., given plausible choices of iuq. Systematic un- 
certainties in the AGN LF do not appear large enough to ac- 
commodate e > 0.15. Accommodating these high efficiencies 
would instead require a substantial downward revision of re- 
cent estimates of the local black hole mass function, reducing 
the integrated mass density to p, 2 x IO^MqMpc"^. Our 
results are consistent with a scenario like the one of King & 
Pringle (2006) in which chaotic accretion spins down black 
holes because of counter-alignment with the accretion disk 
angular momentum, or with other mechanisms that reduce ef- 
ficiencies below the MHD-simulation predictions. 

The assumption that all active black holes at a given mass 
and redshift have the same riio is clearly an idealization, at 
least in the local universe where observations indicate a wide 
range of Eddington ratios (Heckman et al. 2005; Greene & 
Ho 2007). Steed & Weinberg (2003) and Yu & Lu (2004) 
discussed continuity equation models evolved adopting a dis- 
tribution of Eddington ratios. In particular, Yu & Lu (2004) 
have derived the relation between the integrated number of 
AGNs shining at all times at a given luminosity L, the mean 
light curve of black holes, and the local black hole mass func- 
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tion. Following their equation (18), we find that a good match 
between the cumulative number of AGNs and of relic black 
holes can be achieved for e ^ 0.07 (required by the Soltan ar- 
gument) and a mean AGN light curve exponentially increas- 
ing with A = 0.6 and a negligible declining phase, similar to 
our reference model. Alternatively, we find that a good match 
can be obtained by assuming that black holes grow rapidly 
in a Super-Eddington phase with A > 2 and then have a long 
declining p hase, qualitatively resembling our m{z) model dis- 
cussed in § 14.41 In future work we will investigate models 
that incorporate multiple mo-values and accretion modes, in- 
cluding the addition of modest log-normal scatter in mo{M,z) 
(e.g., Kollmeier et al. 2005; Netzer et al. 2007; Shen et al. 
2007b) and sharper revisions in which some black holes ac- 
crete at super-Eddington or highly sub-Eddington rates, per- 
haps with reduced radiative efficiencies (Narayan, Mahade- 
van, & Quataert 1998 and references therein). We will also 
incorporate mergers at the rates predicted by theoretical mod- 
els of cold dark matter subhalos and their associated black 



holes (Yoo et al. 2007). For appropriate parameter choices, 
we expect that many scenarios can be made consistent with 
the observed AGN LF and the local black hole mass function. 
However, direct measurements of Eddington ratio distribu- 
tions and measurements of AGN clustering and host proper- 
ties, all as a function of luminosity and redshift, should greatly 
narrow the field of viable models. Within the (often substan- 
tial) uncertainties of existing data, a simple model in which 
all black holes grow by accreting gas at mildly sub-Eddington 
rates with a radiative efficiency e « 0.06 — 0.1 is surprisingly 
successful at reproducing a wide range of observations. 
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TABLE 1 

AGN Luminosity Function parameters 



Parameter 



Pi 


4.23 (z < 4) 


-0.615Z + 6.690 


3.0 (z = 6) 


P2 


-LSatalU 






71 


0.86 at all z 






72 


2.6 (z < 0.5) 


--0.933Z + 3.067 


2.32 (0.8 < z < 4) 0.24Z+ 1.36 2.8 (z = 6) 


A 


5.04 X 10-'*Mpc-3 






z? 


1.9 






logLn 


44.6 (z < 3) 


0.0667Z+44.4 


44.8 (z = 6) 


logL. 


43.94 







Note. — List of the parameters entering the AGN luminosity function described in §|2] Some of the 
parameters assume different values in different redshift bins as quoted. A smooth linear transition in redshift z 
in the values has been applied between discontinuous redshift bins. Luminosities are the inferred values in the 
2—10 keV band, in units of erg s~ ' . 



TABLE 2 

AGN BOLOMETRIC LUMINOSITY FUNCTIONS 





logL 


log* 


log'J'noCT 


log $DoublcCT 


log *HRH07 


(1) 


(2) 


(3) 


(4) 


(5) 


(6) 


0.020 


41.00 


-1.859 


-1.994 


-1.757 


-1.374 


0.020 


41.25 


-2.074 


-2.208 


-1.972 


-1.578 


0.020 


48.00 


-10.08 


-10.16 


-10.01 


-9.539 


0.260 


41.00 


-1.903 


-2.036 


-1.796 


-1.609 


0.260 


41.25 


-2.089 


-2.222 


-1.984 


-1.787 


0.260 


41.50 


-2.267 


-2.399 


-2.162 


-1.967 



Note. — (1) redshift; (2) bolometric luminosity in logarithmic 
scale in units of ergs"'; (3) reference bolometric AGN luminosity 
function in logarithmic scale in units of Mpc^^ dex^ ' ; (4) reference 
bolometric AGN luminosity function with no Compton-thick sources 
(i. e., with logW/y/cm^^ > 24); (5) reference bolometric AGN lumi- 
nosity function with Double contribution from Compton-thick sources 
(i. e., with logNn /cm~^ < 26); (6) bolometric AGN luminosity func- 
tion from HRH07. The full table is available in electronic form. 



TABLE 3 
BLACK HOLE MASS FUNCTIONS 



z 

(1) 


logM. 

(2) 


log$ 
(3) 


log ^noCT 

(4) 


log 'JDoublcCT 

(5) 


log 'I'HRHO? 

(6) 


log*HRH07EXTRA 

(7) 


0.02 


5.0 


-1.264 


-1.398 


-1.162 


-0.928 


-0.952 


0.02 


5.2 


-1.394 


-1.528 


-1.292 


-1.060 


-1.082 


0.02 


9.6 


-4.683 


-4.779 


-4.604 


-4.436 


-4.964 


0.26 


5.0 


-1.482 


-1.614 


-1.376 


-1.129 


-1.131 


0.26 


5.2 


-1.608 


-1.739 


-1.501 


-1.247 


-1.248 


0.26 


5.4 


-1.723 


-1.855 


-1.616 


-1.365 


-1.364 



Note. — (1) redshift; (2) black hole mass in logarithmic scale in units of Mq ; (3) black 
hole mass function in our reference model in logarithmic scale in units of Mpc~' dex~ ' ; (4) 
reference black hole mass function in the same model but with no Compton-thick sources; 
(5) reference black hole mass function in the same model but with Double contribution from 
Compton-thick sources; (6) black hole mass function in the same model but adopting the 
HRH07 luminosity function; (7) black hole mass function obtained adopting the HRH07 
luminosity function with e = 0.09 and m = 1, as in Figure[8jl. The full table is available in 
electronic form. 
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TABLE 4 
BLACK HOLE DUTY CYCLES 



Z logM, Po Po"°CT pDoubleCT pHRH07 ^RHOTbxtra 



(1) 


(2) 


(3) 


(4) 


(5) 


(6) 


(7) 


0.02 


5.0 


0.00995 


0.00995 


0.00995 


0.01401 


0.00727 


0.02 


5.2 


0.00941 


0.00942 


0.00942 


0.01299 


0.00668 


0.02 


9.6 


0.00008 


0.00008 


0.00008 


0.00011 


0.00009 


0.26 


5.0 


0.03063 


0.03057 


0.03050 


0.01947 


0.01050 


0.26 


5.2 


0.03246 


0.03231 


0.03215 


0.01837 


0.00987 


0.26 


5.4 


0.03273 


0.03243 


0.03212 


0.01729 


0.00924 



Note. — (1) redshift; (2) black hole mass in logarithmic scale in units of 
Mq; (3) duty cycle for our reference model; (4) duty cycle for the reference 
model with no Compton-thick sources; (5) duty cycle for the reference model 
with Double contribution from Compton-thick sources; (6) duty cycle for the 
model obtained adopting the HRH07 luminosity function; (7) duty cycle for the 
model obtained adopting the HRH07 luminosity function with e = 0.09 and rh = 
1. The fuU table is available in electronic form. 
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Fig. 1. — The bolometric AGN luminosity function. Curves sliow tlie model described in §|2] The solid line is the total LF including very Compton- 
thick sources, the short-dashed line includes sources with column densities up to logA'/y/cm"^ = 24, while the dot-dashed line includes only sources with 
logW/y/cm"- < 22. The data from optical surveys are from Hunt et al. (2004), Pei (1995), Wisotzki (1999), Jiang et al. (2006), Cool et al. (2006), Shankar & 
Mathur (2007; derived from data of Stern et al. 2000; Willott et al. 2004; Mahabal et al. 2005), VIMOS-VLT Deep Survey (Bongiomo et al. 2007), 2SLAQ 
(Richards et al. 2005), SDSS (Richards et al. 2006; Fan et al. 2001, 2004), COMBO-17 (Wolf et al. 2003), Kennefick et al. (1994). The data from X-ray surveys 
are from Ueda et al. (2003), Barger et al. (2003), Barger et al. (2005), Barger & Cowie (2005), Nandra et al. (2005) and Silverman et al. (2008). GOODS 
(multi-wavelength) data are from Fontanot et al. (2007). 
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2. — Same as Figure[T]but now all the data have been corrected for the obscured fraction as given in HRH07. 
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log N„<25 •HEAO 




1 10 100 

E (keV) 

Fig. 3. — The X-ray background. Curves represent the integration of our AGN bolometric LF converted to hard X-ray bands, as described in the text. The 
dashed line includes sources with column densities up to logWH/cm^^ < 22, the dot-dashed line those up to logW/y/cm^^ < 23, and the triple-dot dashed 
line those with logW/y/cm^- < 24. The solid line refers to the total when including Compton-thick sources up to logA'/^/cm"^ < 25, with a dependence on 
luminosity following that in U03. Data are taken from a compilation in Frontera et al. (2007). 
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Fig. 4. — Comparison between the bolometric AGN luminosity function described in ^^{solid line) and the recently derived estimate by HRH07 {dashed 
hne), shown with its I-ct uncertainty in the bolometric coiTection (gray area). Curves are vertically offset by —1.5 X z. 
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Log M. [Ms] 

Fig. 5. — Comparison among estimates of the local black hole mass function. Lines show estimates using various calibrations of the M, — Lsph, — or 
M, — Mstai relations as described in the text, assuming a 0.3-dex intrinsic scatter in all cases. The grey band encompasses the range of these estimates. Filled 
small circles show the determination of Hopkins et al. (2()07b) using the black hole "fundamental plane" and open circles show the determination of Graham et 
al. (2007) using the relation between black hole mass and Sersic index. 
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Fig. 6. — Black hole growth and stellar mass growth, {a) Average black hole accretion rate as computed via equation @ compared to the SFR as given by 
Hopkins & Beacom (2006) and Fardal et al. (2007), scaled by the factor M./Mstar = 0.5 x 1.6 x 10"^ The grey area shows the 3-cr uncertainty region from 
Hopkins & Beacom (2006). (b) Cumulative black hole mass density as a function of redshift. The solid line is the prediction based on the bolometric AGN 
luminosity function. The light-gray area is the local value of the black hole mass density with its systematic uncertainty as given in Figure|5] The dark squares 
are estimates of the black hole mass function at z = 1 and z = 2 obtained from the stellar mass function of Caputi et al. (2006) and Fontana et al. (2006), scaled 
by the local ratio M. /Mstar = 1.6 x 10^'. The lines are the integrated stellar mass densities based on the SFR histories in panel (a), scaled by 8 x lO^''. 
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Fig. 7. — Results for a reference model that agrees well with the local black hole mass function, with e = 0.065, nio = 0.60, and an initial duty cycle 
fo(z = 6) = 0.5. (a) The accreted mass function shown at different redshifts as labeled, compared with the local mass function {grey area), {b) Similar to (a) but 
plotting M, $ (M, ) instead of <3? (M, ) . (c) Duty-cycle as a function of black hole mass and redshift as labeled, {d) Average accretion histories for black holes of 
different relic mass M, at z = as a function of redshift; the dashed lines represent the curves when the black hole has a luminosity below Lmin{z) and therefore 
does not grow in mass; the minimum black hole mass corresponding to Lmin(z) at different redshifts is plotted with a dot-dashed line. 
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Fig. 8. — Effect of model inputs on the predicted black hole mass fimction at z = 0. (a) Varying e at fixed njo = 0.6. (2>) Varying = at fixed e = 0.065. (c) 
Solid and dotted curves show the effect of removing or doubling the fraction of Compton-thick sources, keeping e = 0.065 and nio = 0.6 fixed at their reference 
values. Dashed and dot-dashed curves show these modified models with parameter values chosen to reproduce the local mass function, (d) The solid curve 
shows the effect of adopting the HRH07 LF, with e = 0.065 and mo = 0.6. The dashed curve shows a model with the HRH07 LF and parameter values e = 0.09 
andnio = 1. 
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Fig. 9. — Models incorporating our standard luminosity function estimate compared to the properties of the local black hole mass function described in § 14.31 
integrated mass density (p., upper left), peak mass (Mpeak. upper right), width (M, [^2. lower left), and asymmetry (Ap,/p,, lower right). The thick lines 
show model predictions as a function of radiative efficiency e, for several values of rfiQ as labeled. The horizontal gre\> band shows observational estimates of 
the four quantities based on the grey band in Figure[5] Horizontal long-dashed and short-dashed lines show the values derived from the local mass functions of 
Hopkins et al. (2007b) and Graham et al. (2007), respectively. 
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Fig. 10. — Local black hole mass function parameters predicted for different LF inputs. The accretion rate is fixed to mo = 0.75 in all cases. The thick lines 
refer to models with our standard LF estimate (solid line), with no contribution of Compton-thick sources (dotted line), with only Type I AGNs included (dashed 
line), with Type 1 AGNs multiplied by a constant factor of 3 at all redshifts and all luminosities (dot-dashed line), and with the HRH07 LF (long-dashed line). 
Grey band and horizontal thin lines are as in Figure |9] 
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Fig. 1 1 . — Model with e = 0.065 and decreasing mo(z) as given in equation (T7J, compared to our reference model with constant liio = 0.6 and e = 0.065. 
Panels (a), (b), (c) show, respectively, the evolution of the mass function, the evolution of the duty cycle, and mass growth along characteristics, in the same 
format as Figure|2] Reference model results are shown by squares in (a), by thinner lines in (b), and by long-dashed lines in (c). 
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Fig. 12. — Model with e = 0.075 and a mas s-de pendent Eddington ratio mo = 0.445 (M, /lO' Mq)"^ (equation llSt , compared to our reference model. The 
format is similar to panels (a) and (b) of Figure llll 
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Fig. 13. — Evolution of the black hole mass function in a model with black hole mergers. Accretion-driven growth is computed assuming e = 0.065, mg = 0.60, 
and each black hole has a probability Pmcrs = 0.5 of merging with another black hole of equal mass per Hubble time tfi{z)- Squares show the z = predictions 
of the reference model (same accretion parameters, no merging), and the long-dashed line shows the z = mass function for a merger model with the HRH07 LF 
and accretion parameters e = 0.09, »io = 1 ■ 
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Log L [erg/s] Log L [erg/s] 

Fig. 14. — Panel a: cumulative space density of predicted hosts (including those with inactive black holes of the same mass) th at h ave a lum inosity higher than 
L for four different redshifts. Panel b: predicted bias as a function of luminosity and redshift, computed according to equations (22) and )23t : the data are from 
da Angela et al. (2Q06, filled triangles), Myers et al. {2007, filled squares), Porciani & Norberg (2007, filled and open circles), and Shen et al. (2007a; open 
diamonds), all corrected to ug = 0.8. Both panels refer to the predictions of our reference model. 



